home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / scripts / control / are.m < prev    next >
Encoding:
Text File  |  1996-07-16  |  2.7 KB  |  93 lines
  1. ## Copyright (C) 1996 John W. Eaton
  2. ##
  3. ## This file is part of Octave.
  4. ##
  5. ## Octave is free software; you can redistribute it and/or modify it
  6. ## under the terms of the GNU General Public License as published by
  7. ## the Free Software Foundation; either version 2, or (at your option)
  8. ## any later version.
  9. ##
  10. ## Octave is distributed in the hope that it will be useful, but
  11. ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13. ## General Public License for more details.
  14. ##
  15. ## You should have received a copy of the GNU General Public License
  16. ## along with Octave; see the file COPYING.  If not, write to the Free
  17. ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
  18. ## 02111-1307, USA.
  19.  
  20. ## Usage: x = are (a, b, c {,opt})
  21. ##
  22. ## Solves algebraic riccati equation
  23. ##
  24. ##   a' x + x a - x b x + c = 0
  25. ##
  26. ## for identically dimensioned square matrices a, b, c.  If b (c) is not
  27. ## square, then the function attempts to use b * b' (c' * c) instead.
  28. ##
  29. ## Solution method: apply Laub's Schur method (IEEE Trans. Auto. Contr,
  30. ## 1979) to the appropriate Hamiltonian matrix.
  31. ##
  32. ## opt is an option passed to the eigenvalue balancing routine default is "B".
  33. ##
  34. ## See also: balance
  35.  
  36. ## Author: A. S. Hodel <scotte@eng.auburn.edu>
  37. ## Created: August 1993
  38. ## Adapted-By: jwe
  39.  
  40. function x = are (a, b, c, opt)
  41.  
  42.   if (nargin == 3 || nargin == 4)
  43.     if (nargin == 4)
  44.       if (! (strcmp (opt, "N") || strcmp (opt, "P") ...
  45.          || strcmp (opt, "S") || strcmp (opt, "B") ...
  46.          || strcmp (opt, "n") || strcmp (opt, "p") ...
  47.          || strcmp (opt, "s") || strcmp (opt, "b")))
  48.     warning ("are: opt has an invalid value; setting to B");
  49.     opt = "B";
  50.       endif
  51.     else
  52.       opt = "B";
  53.     endif
  54.     if ((n = is_square(a)) == 0)
  55.       error ("are: a is not square");
  56.     endif
  57.  
  58.     if (is_controllable(a,b) == 0)
  59.       warning ("are: a, b are not controllable");
  60.     endif
  61.     if ((m = is_square (b)) == 0)
  62.       b = b * b';
  63.       m = rows (b);
  64.     endif
  65.     if (is_observable (a, c) == 0)
  66.       warning ("are: a,c are not observable");
  67.     endif
  68.     if ((p = is_square (c)) == 0)
  69.       c = c' * c;
  70.       p = rows (c);
  71.     endif
  72.     if (n != m || n != p)
  73.       error ("are: a, b, c not conformably dimensioned.");
  74.     endif
  75.  
  76.     ## Should check for controllability/observability here
  77.     ## use Boley-Golub (Syst. Contr. Letters, 1984) method, not the
  78.     ##
  79.     ##                     n-1
  80.     ## rank ([ B A*B ... A^   *B]) method
  81.  
  82.     [d, h] = balance ([a, -b; -c, -a'], opt);
  83.     [u, s] = schur (h, "A");
  84.     u = d * u;
  85.     n1 = n + 1;
  86.     n2 = 2 * n;
  87.     x = u (n1:n2, 1:n) / u (1:n, 1:n);
  88.   else
  89.     usage ("x = are (a, b, c)");
  90.   endif
  91.  
  92. endfunction
  93.